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Abstract 

We have studied channeling effects in a Cesium Iodide (Csl) crystal that is 
similar in composition to the ones being used in a search for Weakly Inter¬ 
acting Massive Particles (WIMPs) dark matter candidates, and measured its 
energy-dependent quenching factor, the relative scintillation yield for elec¬ 
tron and nuclear recoils. The experimental results are reproduced with a 
GEANT4 simulation that includes a model of the scintillation efficiency as 
a function of electronic stopping power. We present the measured and sim¬ 
ulated quenching factors and the estimated effects of channeling. 
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1. Introduction 

Weakly Interacting Massive Particles (WIMPs) are candidates for dark 
matter [1] that satisfy the relic density of dark matter observed by the 
WMAP |2] and PLANCK [3] experiments and successfully explain the ob¬ 
served gravitational lensing effects of galatic clusters Si- In the Standard 
Dark Matter Halo Model, WIMPs are spread throughout the galaxy with 
a Maxwellian velocity distribution and transfer their kinetic energy to ordi¬ 
nary matter by WIMP-nucleon elastic scattering PE]. The scattered nuclei 
subsequently produce detectable responses in the material via ionization, 
scintillation and phonon creation, from which WIMP-nucleon scattering can 
be inferred. 

Elastic scattering of WIMPs with masses of hundreds or thousands of 
GeV on nuclei would produce recoil nuclei with energies ranging from tens to 
hundreds of keV. Since it is hypothesized to be a weak interaction process, 
it is expected to occur with a very small probability. Thus, a WIMP detec¬ 
tor is required to have a high sensitivity for nuclear recoils and a very low 
level of radioactive background contaminants. Thallium-doped Gsl crystals 
(GsI(Tl)) are well suited for WIMP searches thanks to their high scintilation 
light yields, about 65,000 photons/MeV for electron recoils [SI E], and the 
low radioactive background levels that can be achieved through purihcation 
techniques applied during the growing process mm- A particular advan¬ 
tage of GsI(Tl) detectors is their different time responses for electron- and 
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nuclear-induced ionization, which permits the use of pulse shape discrim¬ 
ination (PSD) techniques to remove signihcant fractions of electron-recoil 
background events induced by radioactive background sources imilS]. Spin- 
independent WIMP-nucleus elastic scattering is a coherent process with a 
cross section that has quadratic dependence on the atomic mass number; 
the cross section for spin-dependent elastic scattering has a quadratic depen¬ 
dence on the nuclear spin expectation value |7]. Because of their large atomic 
mass numbers and nuclear-spin expectation values, Cesium (Cs) and Iodine 
(I) ions are expected to have relatively large WIMP-nucleus cross sections 
for both the spin-independent and spin-dependent scenarios. The Korea In¬ 
visible Mass Search (KIMS) uses CsI(Tl) crystals for WIMP searches and 
has published stringent upper limits for WIMP-proton spin-dependent elas¬ 
tic scattering cross-sections [HE]. 

The quenching factor (QF) is the scintillation light yield produced by 
a nuclear recoil relative to that for an electron recoil at the same energy. 
Typically electron-recoil responses are measured using gamma-ray sources, 
in which case the QF can be expressed as: 


QF 




recoil 


-^meas -S7,calib 

77 r ' 

-^recoil -^7,calib 


( 1 ) 


where F^recoii, E^^as are the recoil energy of the nucleus produced by WIMP- 
nucleus scattering and the experimentally measured energy by the crystal 
detector, respectively. Since the light output, Tmeas, is measured in the scin¬ 
tillator, Ejneas Can be obtained from the calibration. L..y,caiib is the measured 
scintillation light yield for gamma rays with a known energy E^^caiih- 

The quenching factors of Csl crystals used in the KIMS experiment have 
been previously measured na. Subsequently, it has been suggested that 
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channeling effects in the scintillation crystal might enhance the qnenching 
factors to as high as QF~ 1 for some specific nuclear-recoil conditions ng. 
Channeling occurs when a recoil ion in the target material moves in a di¬ 
rection that is within a critical angle from a symmetry axis or plane in the 
crystal lattice [IHlin]. In these cases, the recoil ion primarily loses energy 
via numerous scatterings with atomic electrons around target nuclei that are 
confined to small scattering angles because of the relatively large impact pa¬ 
rameters between the moving ion and target nuclei. As a result, channeling 
effects show up as enhanced ion penetration ranges and larger numbers of 
electron-hole pairs in the target material with a resultant smaller stopping 
power. The large penetration range enhances the scintillation yield in the 
crystal in accordance with Birks’ formula [iHl UHl 1^ . 

However, the aforementioned channeling effects have only been studied 
using ions that are incident from outside of the crystal into the empty space 
between symmetrically aligned lattice atoms na. Thus, conclusions from 
previous studies of channeling effects might not apply to recoil ions from 
WlMP-nucleus scattering, which occur inside the crystal. In these cases, 
recoil ions, originally located at a crystal lattice point, initially travel at a 
large angle from the adjacent target nuclei located near the symmetry axis 
or crystal plane because of a small initial impact parameter. As a result, the 
recoil ions cannot easily channel through the empty space near the symmetric 
axis; this is known as the blocking effect. Bozorgnia et al. [2T1 [ 22 ] point out 
that recoiling lattice ions have some chance to be channeled through the 
symmetric axis due to thermal lattice vibrations. Nevertheless, the fraction 
of full channeling for isotropically scattered recoil ions is expected to be below 
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10 % for Csl crystals [22]. 

In this report, we present new measurements of quenching factors for a 
monocrystalline CsI(Tl) crystal. In addition, we provide estimates of chan¬ 
neling effects in the crystal by comparing the measured energy spectrum 
with the result of the GEANT4 simulations coupled with a program called 
MARLOWE [23] |21|. The impinging neutrons in the experimental setup 
were modelled with a GEANT4-based simulation and the propagation of 
the recoil ions in the monocrystalline structure was incorporated using the 
MARLOWE program. 


2. Experiment 
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Figure 1: Schematics of the experimental setups, (a) Setup I and (b) Setup II. The arrows 
indicate the direction of the incident neutron beam. The numbers next to the neutron 
detectors represent polar scattering angles, 0 nd. The shielding and supporting materials 
surrounding detectors are not drawn for clarity. 


To study nuclear recoil signals in GsI(Tl), mono-energetic 2.4 MeV neu- 
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trons from a neutron generator were irradiated onto a Csl crystal. A model 
MP320 neutron generator from ThermoFisher Scientific produces neutrons 
isotropically from ^H(^H,n)^He reactions. The maximum neutron flux is 
7 X 10^ /s for a 60 /rA beam current of 90 keV deuterons. However, for these 
measurements we used 50 keV deuterons, for which the detected flux of neu¬ 
trons in a test experiment that placed a neutron detector at the output of 
the neutron generator was maximum, while the gamma flux was minimized. 

In order to measure the quenching factor and study channeling effects, 
two different experimental setups were used. The hrst one, called Setup I 
and depicted in Fig. [^(a), was optimized for measuring quenching factors in 
the crystal for various neutron scattering angles, ©nd, where the subscript 
ND indicates the neutron detector. In this setup, neutron detectors were 
located at eight different polar scattering angles, 45°, 60°, 75°, 90°, 105°, 
120°, 135°, and 150°. To increase statistics, at each scattering angle, two 
neutron detectors were placed symmetrically on both sides of the neutron 
beam direction. For 90°, to avoid blocking of the scattered neutrons by the 
PMT attached to the Csl crystal, the neutron detectors were positioned in 
the plane perpendicular to the direction of the neutron. The distances from 
the Csl crystal to the neutron detectors ranged from 15.5 cm to 40.0 cm 
in order to partially compensate for the variation in cross sections of the 
neutron elastic scattering at the different scattering angles. 

The second setup. Setup II, described in detail in ref. |25] and depicted 
in Fig.[g(b), is designed to study channeling effects in nuclear recoil events. 
This has six neutron detectors at a 60° polar angle that cover different az¬ 
imuthal recoil-ion directions. The six neutron detectors are located 1.1 m 
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away from the Csl crystal in order to restrict -Erecoii values to a narrow range. 

The Csl crystal had PMTs (9269QA by Electron Tubes, green-extended) 
attached at each end and was contained in a 1 mm thick copper box. The 
neutron detectors consisted of BC501A liquid scintillator in a 5 mm thick 
double-layered-bottle of teflon and stainless steel read out by a R329-02 PMT 
made by HAMAMATSU. The Csl crystal dimensions were 3 x 3 x 1.4 cm^ 
and the Thallium doping levels were under 0.1 % mole. 

The neutron generator was surrounded by a 10 cm thick polyethylene 
shield with an additional 10 cm thick lead wall on the side facing the Csl 
crystal. Neutrons passing through a 15 cm thick-lead guide from the gener¬ 
ator exited the generator via 3.2 cm diameter apertures in the polyethylene 
and lead walls. In Setup I, to reduce ambient gamma backgrounds, 5 cm 
thick lead blocks were put around the neutron detectors and the Csl crystal 
except along the neutron trajectories. The neutron detectors in Setup II 
were contained in 10 cm-thick polyethylene boxes to shield them from ambi¬ 
ent radioactivity. In this setup, the Csl detector was not shielded along the 
neutron entrance and exit directions. 

Events were triggered by a coincidence between valid signals from the 
Csl crystal and any one of neutron detectors in Setup I that occurs within 
a 2 /is timing window. For Setup II, an additional signal from the neutron 
generator was also required in order to reduce random coincident events. A 
valid signal in the Csl crystal was dehned as at least one photoelectron in 
each PMT within a 2 /rs interval. This was realized with a 400 MHz Fast 
Analog to Digital Converter (FADC400) module made by Notice Korea. The 
thresholds of the neutron detectors were set at the single photoelectron level. 
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Figure 2: (color online) A two-dimensional plot showing the pulse shape discrimination 
in one of the neutron detectors. The x-axis shows the sum of all the collected charges for 
the first 1 /is time window. The y-axis shows the fraction of the charges not contained 
the cluster with the largest summed charge. The black dot points are from the neutron 
scattering experiment and the blue crosses are from gamma ray calibration data using a 
^^Na source. Points above the red solid curve are identified as neutron-induced events; 
points below the curve are rejected as gamma-induced events. 

In the offline analysis, nentron scattering events were selected by exploit¬ 
ing the pnlse shape discrimination (PSD) power of the liqnid scintillator 
as shown in the scatter plot of Fig. where the horizontal axis shows the 
snmmed charge over 1 /rs and the vertical axis shows the fraction of the 
charge that is not in the dominant peak, which we call the cluster, with the 
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largest summed charge. The events above (below) the red solid line in Fig. 
are primarily due to neutrons (gamma rays). In the hgure, the small dots 
are data taken with the neutron generator on and the blue crosses are the 
data taken with a ^^Na source located between the Csl crystal and neutron 
detectors. Two 511 keV gamma rays simultaneously emitted from the ^^Na 
source are detected by the Csl crystal and one of the neutron detectors. From 
these data we estimate that only about 0.3 ± 0.1(stat.)% of gamma-induced 
events are misidentihed as neutron events. 

Another criteria for the event selection, called £t quality, required that the 
signal of the Csl exhibited a decay time that is characteristic of Csl crystals 
when fitted with a single exponential function. This criteria removed low 
energy events triggered by the long tails of a high-energy signal in the Csl 
crystal that happened to occur in accidental coincidence with a background 
signal in the neutron detector. 

The temperature of the experimental environment was controlled to be 
stable at 18 °C and an energy calibration with an source was performed 

every three days during the data taking. 

3. Simulation 

A full simulation of the experimental setup was carried out and the re¬ 
sults were compared with the experiment data. The simulation consists of 
two parts: the scattering of the neutron beam in the Csl crystal and the 
subsequent detection in the neutron detector that is simulated by GEANT4. 
At each step, the deposited energy is converted to a scintillation light output 
in keVee units based on a simulation study for recoil ions in a Csl crystal us- 
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(b) Stopping power for an event 


(a) Average stopping power 




Figure 3: The electronic stopping power distribution versus the penetration depth. These 
histograms were obtained from the averaged ionization energy loss at each penetration 
depth bin of the primary I ion and its recoil ions for all simulated events (a) and for a 
typical single event (b). The kinetic energy of the primary I ion was initially set as indicated 
in the figure legends : mono, and amor, indicate that a monocrystalline CsI(Tl) or an 
amorphous CsI(Tl) crystal was used as the target materials in the simulation, respectively. 
The gray line histogram in (a) was constructed with a SRIM simulation; the others were 
obtained with a MARLOWE simulation. 


ing the MARLOWE simulation code with a modihed Birks’ formula [2S]. A 
Monte-Carlo (MC) method is used to construct photoelectron signals for each 
event that is tuned to reproduce the measured energy spectra with an energy 
resolution and a detection efficiency that were compatible with experimental 
observations. 

MARLOWE is a binary cascade simulation code for atomic collisions [23| 
123]. Even though it is not as widely used as SRIM |26] it has the feature of 
simulating ion transport not only in monocrystalline or polycrystalline ma¬ 
terials, but also in amorphous materials. The MARLOWE program reads 
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an input file that sets parameters that are relevant to the crystal structure, 
material properties, physics models for nuclear and electronic stopping, se¬ 
lection criteria for target atoms, etc. Using parameter values taken from the 
papers cited in ref. [2S] with some adjustment, we hnd that the ions range 
and stopping power obtained from the MARLOWE program are compatible 
to those determined from SRIM for amorphous CsI(Tl). The only adjust¬ 
ment we made was to the fraction of the nonlocal part in the Oen-Robinson 
electronic stopping power model [23] • In this model, the nonlocal part is 
independent of the impact parameter between a projectile ion and a target 
atom, while the local part has dependence on the electron density of the 
target atom. The nonlocal fraction can be different for various materials. In 
ref. |27|, the authors used 0.4 to model the ranges of boron ions in silicon. 
In this work we use 0.65. 

Figure shows the electronic stopping power at each penetration depth 
for Iodine (I) ions obtained from the averaged ionization energy loss in each 
bin. Each individual event has large fluctuations as shown in Fig. |^(b), but 
the average of all events is the smooth curve shown in Fig. [^(a). For 50 keV 
Iodine ions passing through an amorphous crystal, the average electronic 
stopping power and total penetration range from the MARLOWE simula¬ 
tion are similar to those from SRIM with the parameters mentioned above. 
However, when the target is changed to a monocrystalline structure with 
the target atoms moving with thermal vibrations that are centered at each 
lattice site, the shape of the averaged electronic stopping power distribution 
becomes more asymmetric and the maximum penetration range increases as 
shown in Fig. [^(a). 
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For light recoil ions, such as those commonly considered in scintillation 
efficiency studies [2H], the stopping power is dominated by the electronic 
stopping power. However, for heavy ions, such as Cs or I in a Csl crystal, 
with recoil energy below some level for each nucleus, the nuclear stopping 
power due to phonon creation at each penetration depth bin is larger than the 
electronic stopping power. Since the scintillation light output is generated 
from the electronic stopping power, we used the electronic stopping power 
given by the MARLOWE program as an estimate of the scintillation yield. 

Once the electronic stopping power is obtained from the MARLOWE 
program, we can estimate the differential scintillation light yield per energy 
deposit from Birks’ formula [IS] : 


dL _ S 
dE ~ 1 + 

dr 


( 2 ) 


where L is the generated scintillation yield, E is the energy deposited by 
the ion in the material, S is the absolute scintillation yield per unit deposit 
energy, kB is the Birks’ factor, and ^ is the ionization energy per unit 
penetration depth. For stopping powers below 20 MeV-cm^/g, we used a 
modihed Birks’ formula, based on Eq. (11.8a) in ref. [18], which hts the data 
better [28], 1^ : 

dL ^ ka§ S 

dE 1 + ka^l + kB^’ ' 

where the additional term containing ka describes the correlation between 
the scintillation yield and the number of electron-hole pairs in the low stop¬ 
ping power region. 

As discussed in ref. [25], the parameters in the modihed Birks’ formula 
were obtained by htting Eq. ([^ to the measured scintillation efficiency data [22] 
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for the stopping powers of light charged particles in CsI(Tl) crystals. The 
parameters determined from this £t are ka = 1 g/MeV-cm^, kB = 3.8 x 
10“^ g/MeV-cm^ and S = 1.375 keVee/keV. Equation ^ was then multi¬ 
plied by the average ionization energy at each penetration depth and summed 
for the whole bin to obtain the light output for the Csl crystal as shown in 


Eq. Q. 


max.bin max.bin 



L = E = E 


2 = 0 2=0 


where L is the total light output, and AL* and AEi are the light output and 
the ionization energy at each penetration depth bin. 

The data used to £x the parameters of the modihed Birks’ formula were 
taken with an experimental setup different from the current one. The light 
output produced by gammas depends upon the DAQ time window because 
of the presence of slow component of the scintillation, and this leads conse¬ 
quently to different energy resolutions and some non-linearity I3D1. On the 
other hand, the scintillation from alphas has smaller slow components than 
that for gammas and is, therefore, less dependent on the DAQ time window. 
As a result, the measured quenching factor for alphas can depend on the 
gamma reference line that is used in an experiment. The simulated output 
of Eq. ([^ has been adjusted by an additional correction factor of 0.93 in order 
to reproduce the alpha quenching factors that are consistent with measure¬ 
ment [25]. For the current experiment, the correction factor may be different 
from the one we used, so we introduced another factor, called the shift factor, 
while retaining the scintillation efficiency function multiplied by 0.93. This 
shift factor will be applied to the number of photons in the Monte-Carlo 
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method explained later to vary the relative light outputs of nuclear recoils. 


(a) (b) 




Figure 4: (color online) Distributions of simulated scintillation light output i?meas from 
the simulation of Cs ions with 40 keV recoil energy, (a) is for the monocrystalline CsI(Tl) 
crystal and (b) is for the amorphous Csl. The solid points are the result of the simu¬ 
lation, the red solid lines are the htted results for the Landau-Gaussian functions, and 
the blue dashed lines are the Landau distributions constituting the htting functions. The 
corresponding ht parameters are shown in the inset for each panel. 


The solid poiuts iu Fig. show the simulated euergy distributiou (F^meas) 
for Cs ious with a recoil euergy of 40 keV iu mouocrystalliue (a) aud amor¬ 
phous CsI(Tl) material (b) usiug the MARLOWE program with the modi- 
hed Birks’ formula. Here we varied the iuitial directious of ious isotropically 
withiu ±20° with respect to the directiou perpeudicular to the target plaue. 
The iuitial positious of the Cs ious were set at lattice sites to simulate the 
motiou of recoil ious iuside the crystal. We ht the simulatiou results usiug 

^Erratum : DAQ time window of 25 /iS should be changed to 10 /iS on page 3 in ref. |25) . 
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Landau-Gaussian functions for the Probability Density Functions (PDFs): 

PDF = —Landau{E - El) ^ / 2 a}.)^ 

Cl y 21X00, 

where El is the most probable value of the Landau distribution, Eq is the 
mean value of the Gaussian distribution for the measured energy, and o\, and 
og are the standard deviations of the two distributions. This convolution 
function is calculated by the ROOFIT module in the ROOT program [31] 
and Landau is the Landau distribution function in the ROOT math library. 
The result of these hts are shown as solid red curves in Fig. 

As is evident in Fig. the shape of the PDF for the monocrystalline is 
different from that for the amorphous crystal. The enhancement of events 
in high energy tail region for the monocrystalline case were correlated with 
the ions’ ranges, so these events are identified as being due to channeling 
events. These are from recoil ions that are captured in a channel with the 
help of lattice vibrations and/or the reduction of their energies after several 
scatterings. The lower the ion energy is, the larger the critical angle or the 
probability of the channeling. According to Hobler [32], this continues until 
the critical approach distance between a projectile ion and a target atom is 
similar to the channel radius, the half distance between symmetric axes or 
planes, after which the critical angle goes to zero. Below this energy, ions 
with any incident angle scatter out of the channel, i.e., dechannel. These 
correspond to partially channeled events. Full channeling occurs when ions 
enter a channel at the start of motion and do not escape from the channel 
until they stop. Since an amorphous crystal has no symmetry axes or planes, 
neither partial nor full channeling occurs. 

As mentioned above, we used GEANT4 to simulate the crystal response 
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to neutron-recoil events and the quenching and channeling effect were sim¬ 
ulated by the MARLOWE program with the modified Birks’ formula. The 
simulation incorporated the experimental geometries of Fig. [T] Since the 
application of the MARLOWE simulation at each GEANT4 step inside the 
Csl crystal is CPU intensive, we accelerated the simulation by producing a 
template of PDFs as a function of (Erecoii,Emeas) beforehand. We performed 
the full MARLOWE simulation for a series of recoil energies relevant to our 
experimental setup and determined the parameters of the PDF of Eq. ([^ 
for each recoil energy. Then, by fitting these parameters as a function of the 
recoil energy, a template PDF as a function of Erecoii was generated. For each 
neutron-ion scattering in the GEANT4 simulation, the PDF is applied for the 
corresponding recoil energy, and the Fmeas values is chosen randomly with 
probabilities provided by PDF. In the simulation, we used GEANT version 
4.9.6.p02 with the NeutronHP model and G4EmLivermorePhysics builder for 
the nuclear- and electron-recoil simulations [33] • We excluded events that do 
not conserve energy-momentum in neutron-target atom inelastic scattering 
events, and removed some elastic events with some reduction in efficiency. 
Figurej^shows recoil energy spectra deposited in a GsI(Tl) crystal from single 
scatterings of neutrons on target atoms obtained from GEANT4 simulation, 
corresponding to neutrons that, after the scattering, enter one of neutron 
detectors and deposit some energy. Solid lines depict recoil energy spectra 
from single elastic scatterings and dashed lines are those from single inelastic 
scatterings. Most of the inelastic scattering events in these spectra corre¬ 
spond to those in which gamma rays escape. However, high energy peaks 
from additional energy deposited by gamma rays from excited nuclei can be 
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seen in several figures. For the energy spectra from single elastic scatterings, 
there are unexpected excesses of entries around the expected recoil energy 
peaks. These events are generated by multiple scatters in the shielding ma¬ 
terials after a single scattering in the CsI(Tl) crystal. The relative fraction 
of multiply scattered events increases as scattering angles where the cross- 
section of neutron-nucleus single elastic scattering is small. Because of this, 
we select the single hit events in the GEANT4 simulation for the estimation 
of quenching factors by limiting the recoil energies to a Region of Interest 

(ROI) • (E/j-ecoil,inelastic) 3(7j-ecoil,inelastic ^ E/fecoil ^ (E/j-ecoil,elastic) T SUrecoil,elastic; 

where (Erecou) and (Jrecoii are the Gaussian mean and sigma of single inelastic 
or elastic scattering peak, respectively. 

After we determined the measured energy, or the electron equivalent en¬ 
ergy, from a GEANT4 simulation, we applied a Monte-Garlo method to gen¬ 
erate photons corresponding to that energy. In this Monte-Garlo process, 
we chose a random value that follows a Poisson distribution for a photon 
yield whose mean number per keV is determined by the gamma calibration. 
For energy deposited by gamma rays, the resolution was equivalent to that 
of gamma calibration data using an ^^^Am source, while for recoil ions, the 
resolution was a combined result of a random choice from a Poisson distri¬ 
bution and the photo-electron charge distribution in the experiment. After 
the Monte-Garlo data are generated, we apply the two-fold trigger condition 
requirements as a hardware cut and the £t quality requirements as a soft¬ 
ware cut to mimic those that were applied to the experimental data. This 
Monte-Garlo process is repeated with the application of a shift factor to the 
photon yield from a recoil ion determined from the modified Birks’ formula 
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and the gamma calibration. Since the htted formnla is based on the other 
experimental data, the repetition is a way to hnd a fnnction that describes 
onr data better. In the next section, we present the best £t spectra 

for data and simnlation at each setnp of the scattering angle nsing the most 
probable shift factor obtained from the minimnm negative log likelihood £t 
method. 

In snmmary, we have simulated the distributions of measured CsI(Tl) en¬ 
ergy for several recoil energies using a MARLOWE-based simulation together 
with a modified Birks’ formula Eq. @ [HIES]. The resultant distributions 
are fitted with Landau-Gaussian functions that are incorporated into the 
GEANT4 simulation. And then, by a Monte-Garlo method the signal shape 
per each event was reproduced, and the simulated Emeas spectra were com¬ 
pletely reproduced after applying the selection-efficiency corrections. 

This method to obtain Emeas for each Erecoii can be compared to the work 
of Bernabei et ah [12], where from channelling events was obtained 

using a simple model without full simulation. Previous work by Hitachi [T9| 
and Tretyak |20| simply used the stopping power and the scintillation effi¬ 
ciency function to obtain the mean of E^^as or, equivalently, the quenching 
factors. The work reported here differs from the previous work in that the 
simulation of neutron-ion scatterings and ion-ion scatterings are fully incor¬ 
porated using the GEANT4 and MARLOWE programs, so the simulated 
distributions can be directly compared with experimental measurements of 
both elastic and inelastic scattering, including channeling effects. 
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4. Result and discussion 


Figure]^ shows the experimentally measured (shaded histogram) and the 
simulated (dashed line) -Fmeas distributions for nuclear recoil events tagged 
by each neutron detector in experiment with Setup I. 

Low energy peaks near 1 keV in the simulated spectra in Fig. [^(e)-(h) 
are found to come from neutrons that scatter in the lead blocks around the 
Csl detector after scattering in the crystal. For Fig. |^(a)-(d), those events 
contaminate the single hit events with recoil energy in the ROI, which is 
reflected in discrepancies between the red solid curves denoting the simulated 
spectra of single hit events within ROI and the dashed lines that indicate the 
distribution of all of the simulated events. Furthermore, bumps in right 
hand side of the peaks from single elastic scatterings shown in Fig. [^(e) 
and (f), the multiply scattered events with shielding material, distort the 
measured energy distributions, since their i^meas are within the sigma values 
of the Gaussian distributions for the single elastic events. This implies that 
the quenching factors estimated from the F^meas spectra from the experiment 
may contain systematic errors. So in order to estimate the correct quenching 
factor from the simulation, we consider only single hit MC events within the 
ROI, which are shown as solid red histograms in Fig. Since we applied a 
stronger condition to select the single hit events for the setup of the ©nd = 
90° due to larger solid angle of incident neutrons than in other experimental 
setups, the fraction of single hit events within the ROI is smaller than for 
other cases in Fig. [^(d). 

As mentioned in the previous section, the most probable shift factor for 
the photon yield was applied in order to correctly reproduce the F^meas fitted 
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value from the data. For eight different ©nd values, the mean and sigma of 
the most probable shift factor was 1.133 and 0.105. For ©nd = 75° and 105° 
in Fig. [^(c ) and (e), the shift factors are larger that the rms spread (a). 
These fluctuations may be due to the contamination of multiply scattered 
neutrons in the lead blocks around neutron detectors that are not included in 
our GEANT4 simulation, and the relatively low cross-sections of the neutron- 
nucleus elastic scattering at those angles. Actually we expect the shift factor 
might be smaller than the one we used because of non-linearity between the 
59.54 keV reference line we used, and 662 keV gammas used in ref. 
However, there may be additional factors, for example, differences in the 
relative scintillation efficiencies for alpha and nuclei [2S]; these will be studied 
later. 

The black circles in Fig. [^(a) depict the measured raw quenching fac¬ 
tors from this experiment. For the calculation of the quenching factors, we 
set the denominator in Eq. (|^ to the average Erecoii within the ROI dis¬ 
cussed in section [3] with the error taken to be its standard deviation. The 
numerator is the Gaussian mean value of Emeas in the data, where the £t 
range was determined by the category of the normal scattering represented 
in Table [Tj The quenching factor from MARLOWE was represented as a 
spline function using the Gaussian £t mean for the Emeas spectrum shown in 
Fig. §( a) and given by the MARLOWE program with the modihed Birks’ 
formula for mono-energetic ions. This shows energy dependent quenching 
factor without considering the shift factor in our scintillation model. With 
the red histograms in Fig. [^from the selected events within ROI in the repro¬ 
duced Emeas spectra after applying the shift factor, we obtained the corrected 
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quenching factors as red square points in Fig. The lower value at ©nd = 
75° and higher values at ©nd = 105° are conjectured to be due to a larger 
contamination of multiply scattered neutron events in the lead blocks around 
the neutron detectors that is not included in our GEANT4 simulation, and 
associated with an initial scattering angle in the Csl crystal that is different 
from the nominal value. 

Figure shows the distribution of measured energies, Emeas, obtained in 
Setup II with the same Csl crystal. Here the contributions from the simulated 
-Emeas due to single elastic scatters, multiple scatters, and inelastic scatters in 
the Csl crystal are indicated separately. The single elastic scattering events 
are dominant in the low energy peak region above 1 keV. From the hgure, 
it is clear that the simulation result with the amorphous Csl crystal cannot 
explain the experimental data (Fig. [^(c)), while the simulation with the 
monocrystalline Csl reproduces the experimental data well for all energy 
region (Fig. |^(a)) as expected in Fig. 

Figures [^(b),(d) show the measured energy distributions over an ex¬ 
tended energy range, where two inelastic scattering peaks; 57.6 keV gamma 
rays from the deexcitation of the hrst excited state of I and 79.6 and 81 keV 
gammas from the deexcitation of the hrst excited state of Cs, with added con¬ 
tributions of the recoil energy, can be seen. These experimental distributions 
are compared with the simulation after normalizing by the total number of 
events in a limited energy region of normal scattering as indicated in Table 
and used in Fig. There is an excess of events in Fig. |^(b),(d) above 30 keV 
in x-axis that we attribute to gamma-induced events in the neutron detector 
that are misidentihed as neutron events by the PSD; when we consider the 
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Table 1: Event categories depending on the Emeus 


Gategory 

Energy range 

Normal scattering 

-^meas -^peak “1“ ScTpeak 

Partial channeling 

-^peak “1“ SCTpeak ^ -^meas ^ (-^recoil) 3(Ti-ecoil 

Full channeling 

(-^recoil) ScTj-ecoil -^meas (-^recoil) “1“ ScTrecoil 


neutron PSD data shown in Fig. the total number of gamma ray events 
in data was about 15,000 per each neutron detector and with consideration 
of the 0.3 ± 0.1(stat.)% PSD misidentihcation fraction determined from the 
^^Na calibration run, the origin of excess of events in F^meas can be accounted 
for. The insets in Fig. [^(b),(d) show events in the range from 10 keV to 
30 keV around the full channeling region, where the simulation with an as¬ 
sumption that the detector crystal is amorphous hardly reproduce data due 
to the dehcit of channeling events in the simulation. 

Table shows the event categories that we used for the estimation of 
quenching factors and channeling effects and they are indicated by arrows 
below the x-axis in Fig. [^(a). cxpeak and cxrecoii are the standard deviations of 
the peak of normal nuclear recoils in the experimental data and the peak of 
the F^recoii In the GEANT4 simulation, respectively. Since the nuclear recoil 
peak cannot be htted successfully with a simple Gaussian ansatz due to the 
asymmetric shape of the data, we set the mean value of E-mens to be Epeak- 
We estimated the partial and full channeling fractions with statistical errors 
in the data shown in Fig. |^as 17.04 ± 1.07(stat.)% and 0.75 ± 0.20(stat.)% 
of the total events for all three categories. However, there are contamina- 
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tions from gamma-induced events, multiply scattered events and inelastic 
scattering events as the simulated spectra shown in Fig. The partial and 
full channeling fractions estimated from the simulation are 5.26% and 0.55% 
after removing multiple and inelastic scattering events in the Csl crystal and 
single hit events with multiple scatters in the shielding material that are es¬ 
timated from the energy spectrum for the amorphous Csl crystal simulation 
shown in Fig. |^(c). In any case, these numbers are too small to affect the 
quenching factor. 

In the simulation, we used two tunable parameters. One is the nonlocal 
content for matching electronic stopping power distributions for each pen¬ 
etration depth for SRIM and MARLOWE. And the other is the correction 
factor for matching measured energy spectra of measurement and simula¬ 
tion, which is resultantly 5% increase in the scintillation efficiency function 
from the calculation that 0.93 is multiplied by the mean shift factor of 1.133. 
Although there are uncertainties in this model, it is a good way to estimate 
the channeling effect on the quenching factor by means of simulation for the 
motions of recoil ions in a crystal. 

5. Conclusion 

We measured the quenching factors of a CsI(Tl) crystal to study channel¬ 
ing effects and to analyse the recoil energy spectra expected from WIMP in¬ 
teractions. By comparing experimentally measured Emeas distributions with 
simulations, we find that the simulation method of Emeas using a scintilla¬ 
tion efficiency function of the electronic stopping power in the Csl crystal 
reproduces the data. The quenching factors reported here are lower than the 
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previous measurements. The discrepancies could be due to the use of differ¬ 
ent gamma sources for the energy calibration and the trigger inefficiencies 
for nuclear recoils at low energy |37] due to smaller photon yield in the pre¬ 
vious measurements. More studies on the quenching factors are in progress 
that address these questions. We were able to estimate the fraction of par¬ 
tial channeling as 5.26% with the simulation, contributing to a tail of events 
with enhanced i^meas, while the 0.55% fraction of fully channeled events is at 
a nearly insignihcant level. While this value may have a large uncertainty, 
however consistent with the result of Bozorgnia et ah that the maximum 
fraction for the full channeling is about 2% at a recoil energy of 20 keV at 
room temperature. 
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Figure 5: Recoil energy spectra for single elastic (solid lines) and single inelastic (dashed 
lines) scattering events of neutrons on Cs or I ions. After fitting each peak with a Gaussian 
function, we determined the interest region recoil energy according to the criteria men¬ 
tioned in the text, and selected events in the measured energy spectra in the simulation 
to calculate the quenching factor. Only those events have the correct correlation between 
the angle of the detected neutron and the energy deposited in the Csl crystal by a nuclear 
recoil. 
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Figure 6: (color online) i?meas distributions for nuclear recoil events tagged by neutron 
detectors of different neutron scattering angles of Setup I (shaded histogram) and the 
simulated distributions from the GEANT4 simulation (dashed lines) below 12 keV. The 
red histograms show Emeas for single MC events that are within the recoil energy regions 
of interest (ROI). 

























































































































(a) Quenching Factor from this work (b) Comparison with other experiments 
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Figure 7: (color online) The CsI(Tl) quenching factors as a function of recoil energy, (a) 
shows the quenching factors of this work, and (b) shows the comparison of quenching 
factors of this work and previous experiments. The black circles and red squares in (a) 
are obtained from fitting Gaussian functions to the experimental i?meas distributions and 
those from red solid lines in Fig. The red squares are shifted by -2 keV along the x-axis 
for clarity. The line in (a) is a spline function obtained from simulation quenching factor 
from MARLOWE simulation depicted in Fig. 
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(a) Assuming Mono. CsI(Tl): 0-22 keV 


(b) Assuming Mono. CsI(Tl) : 0-100 keV 
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(c) Assuming Amor. CsI(Tl): 0-22 keV (d) Assuming Amor. CsI(Tl): 0-100 keV 



Figure 8: (color online) Comparison between the experimental and simulated Ameas dis¬ 
tributions for nuclear recoil events tagged by neutron detectors at 0nd = 60° (shaded 
histograms) for Setup II. The simulated distributions (black open histograms) are decom¬ 
posed into single elastic (narrow blue line), inelastic (dark-gray shaded), multiple scatter¬ 
ings (light-blue dashed line) and unidentified events (green dash-dotted line). (a) is for the 
monocrystalline CsI(Tl), (c) is for the amorphous CsI(Tl), and arrows below the x-axis in 
(a) indicate approximate energy regions categorizing events according to Table(b) and 
(d) are zoomed views of the a;-axes of (a) and (c) and insets in (b) and (d) show events 
around the full channeling region. 
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